Damage Spreading and Criticality in Finite Random Dynamical Networks 
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We systematically study and compare damage spreading at the sparse percolation (SP) limit for 
random boolean and threshold networks with perturbations that are independent of the network size 
N. This limit is relevant to information and damage propagation in many technological and natural 
networks. Using finite size scaling, we identify a new characteristic connectivity K s> at which the 
average number of damaged nodes d, after a large number of dynamical updates, is independent 
of N. Based on marginal damage spreading, we determine the critical connectivity K" parse {N) 
for finite N at the SP limit and show that it systematically deviates from K c , established by the 
annealed approximation, even for large system sizes. Our findings can potentially explain the results 
recently obtained for gene regulatory networks and have important implications for the evolution of 
dynamical networks that solve specific computational or functional tasks. 
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Random boolean networks (RBN) were originally in- 
troduced as simplified models of gene regulation 0, 
focusing on a system-wide perspective rather than on 
the often unknown details of regulatory interactions 3 . 
In the thermodynamic limit, these disordered dynamical 
systems exhibit a dynamical order-disorder transition at 
a sparse critical connectivity K c Q ; similar observations 
were made for sparselyconnected random threshold (neu- 
ral) networks (RTN) BSB- 

For a finite system size N, 
the dynamics of both systems converge to periodic at- 
tractors after a finite number of updates. At K c , the 
phase space structure in terms of attractor periods Q, 
the number of different attractors § and the distribu- 
tion of basins of attraction is complex, showing many 
properties reminiscent of biological networks [2| . 

Often, one is interested in the response of dynamical 
networks to external perturbations; because these sig- 
nals can disrupt the generic dynamical state (fixed point 
or periodic attractor) of the network, they are usually 
referred to as "damage." This type of study has numer- 
ous applications, e.g., the spreading of disease through a 
population II lL 12j | , the spreading of a computer virus on 
the internet [13|, failure propagation in power grids 14j, 
or the perturbation of gene expression patterns in a cell 
due to mutations [HI- Mean-field approaches, e.g., the 
annealed approximation (AA) introduced by Derrida and 
Pomeau [J| , allow for an analytical treatment of damage 
spreading and exact determination of the critical connec- 
tivity K c under various constraints 16|, 17 1. It has been 
shown that local, mean-field-like rewiring rules coupled 
to order parameters of the dynamics can drive both RBN 
and RTN to self-organized criticality [H, 0, E3] ■ 

Mean-field approximations of RBN/RTN dynamics 
rely on the assumption that N — > oo and study the 
rescaled damage d(t)/N (where d(t) is the number of 
damaged nodes at time t). For an application to real- 



world problems, these limits are often not very relevant. 
Going beyond the framework of AA, a number of re- 
cent studies focus on the finite-size scaling of (un-)frozen 
and/or relevant nodes in RBN with respect to N [2ll. |22T| ; 
only few studies, however, consider finite-size scaling of 
damage spreading in RBN [l5|,[23j]. Here, of particular 
interest is the "sparse percolation (SP) limit" [23j, where 
the initial perturbation size d(0) does not scale up with 
network size N, i.e., the relative size of perturbations 
tends to zero for large N. This limit applies to many 
of the above-mentioned real- world networks (e.g., the 
spread of a new computer virus on the internet launched 
from a single computer). In this letter, we systemati- 
cally study finite-size scaling of damage spreading in the 
SP limit for both RBN and RTN. We identify a new 
characteristic point K s , where the expectation value of 
the number of damaged nodes after large number of dy- 
namical updates is independent of N . By the definition 
of marginal damage spreading, we introduce a new ap- 
proach to estimate the critical connectivity K C {N) for 
finite N, and present evidence that, even in the large 
TV limit, the critical connectivity for SP systematically 
deviates from the predictions of mean-field theory. 

First, let us define the dynamics of RBN and RTN. 
A RBN is a discrete dynamical system composed of N 
automata. Each automaton is a Boolean variable with 
two possible states: {0, 1}, and the dynamics is such that 



F:{0,ir~{0,l} 



i N 



(1) 



where F = (/i, fjsr), and each fi is represented 
by a look-up table of Ki inputs randomly chosen from 
the set of N automata. Initially, Ki neighbors and a 
look-table are assigned to each automaton at random. 

An automaton state o\ G {0,1} is updated using its 



2 



corresponding Boolean function: 

a t+1 = f (x l x l x l ) 



(2) 



We randomly initialize the states of the automata (initial 
condition of the RBN) . The automata are updated syn- 
chronously using their corresponding Boolean functions. 
The second type of discrete dynamical system we study 
is RTN. An RTN consists of N randomly interconnected 
binary sites (spins) with states o~i = ±1. For each site i, 
its state at time i+1 is a function of the inputs it receives 
from other spins at time t: 



CT i (t + l)=sgn(/ i (t)) 



with 



N 



fi(t) =^2djO-j(t) + h. 



(3) 



(4) 



The N network sites are updated synchronously. In the 
following discussion the threshold parameter h is set to 
zero. The interaction weights take discrete values 
Cij = +1 or —1 with equal probability. If i does not 



receive signals from j, one has Cj 



0. 




FIG. 1: Average Hamming distance (damage) d after 200 sys- 
tem updates, averaged over 10000 randomly generated net- 
works for each value of A, with 100 different random initial 
conditions and one-bit perturbed neighbor configurations for 
each network. For both RBN and RTN, all curves for different 
N approximately intersect in a characteristic point K s . 

Results. We first study the expectation value d of dam- 
age, quantified by the Hamming distance of two different 
system configurations, after a large number T of system 
updates. Let AT be a randomly sampled set (ensemble) 
of zn networks with average degree K, T n a set of zr ran- 
dom initial conditions tested on network n, and I' n a set 
of zi random initial conditions differing in one randomly 
chosen bit from these initial conditions. Then we have 
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FIG. 2: Upper panels: d as a function of TV for different A: 
7? = 1.0 (+), A' = 1.5 (x), A = 1.8 (RBN, *) and A = 1.7 
(RTN, *), A = 2.1 (□, RBN) and A = 1.9 (□, RTN), A = 2.4 
(A, RTN) and A = 2.6 (o). The lines are fits of Eq. ©to the 
data. Lower panel: Scaling exponents 7(A) as a function of 
A, as obtained from fits of Eq. ©for RBN (+) and RTN (x). 
The dashed/dotted lines mark the asymptotes as discussed in 
the text. 



where df(T) is the measured Hamming distance after T 
system updates. Fig. [1] shows d as a function of the 
average connectivity K for different network sizes N by 
using a random ensemble for statistics. For both RBN 
and RTN, the observed functional behavior strongly sug- 
gests that the curves approximately intersect at a com- 
mon point (K s ,d s ), where the observed Hamming dis- 
tance for large t is independent of the system size TV. 

To verify this finding, let us now study the finite size 
scaling behavior of d in this (SP) limit. For K — ► and 
for large A", it is straightforward to estimate the asymp- 
totic scaling. In the case K — > 0, non-zero damage can 
only emerge if the initial perturbation hits a short loop 
of oscillating nodes (most likely a self-connection or a 
loop of length two, longer loops can be neglected). The a 
priori probability to generate these loops is ~ \/N 2 , and 
their number is proportional to the total number of links, 
KN. Hence, we expect d ~ KN/N 2 oc 1/JV. For large 
K, damage percolates through the system, consequently, 
avalanche sizes are bounded only by the size of the sys- 
tem, and we expect d ~ N. At criticality, the frozen 
core of the network always remains undamaged; for this 
reason, d should be limited by the number of unfrozen 
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FIG. 3: Upper panels: d(K) in semi-log scale, as obtained 
from high- precision simulations near K c (ensemble size: 50000 
random networks with 100 random initial conditions for each 
data-point, transient time: 5000 updates). Lines are fits of 
Eq. [8] Lower panels: Scaling exponents 7 derived from equat- 
ing Eq. [S]and|5] with corresponding errorbars. The intersec- 
tion with 7 = (dashed line) defines K s . 
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FIG. 4: The critical connectivity K° parse (N) in the SP limit 
as a function of TV, calculated from Eq. 1121 Curves are 
power-law fits according to Eq. 1131 straight dashed lines mark 
K-anneaivd and Kg for RBN and j^-pN, respectively. 



nodes n u . Recent studies predict n u ~ TV 2 / 3 at K c for 
Boolean networks [2l|, hence, we expect d < TV 2 / 3 at K c . 
For arbitrary K, we make the scaling ansatz 



d(K, N) = a{K) ■ + d (K), -1 < 7 < 1. (6) 

The two upper panels of Fig. illustrate that for both 
RBN and RTN, the numerically measured finite size scal- 
ing for < K < 3 obeys this scaling ansatz very well. In 
both cases, at some K s slightly below K c , we find a tran- 
sition of 7 from negative to positive values (Fig. [2J lower 
panel). The exact determination of the point (K s ,d s ) 
where 7 « is difficult; because of the slow emergence 



of large damage events near K c , measurements with fi- 
nite T can substantially underestimate d (in particular 
for N > 512). We performed high precision numeri- 
cal experiments in the interval 1.6 < K < 2.1, waiting 
for T — 5000 update time steps to let the network dy- 
namics relax after the initial one-bit perturbation; these 
simulations conclusively show an exponential dependence 
d cx exp(c(A) • K) in this interval, with a constant c(N) 
depending only on N (Fig. [3l upper panels). This expo- 
nential dependence becomes apparent with the following 
assumptions: an increase Ad of the average damage is 
proportional to d itself (damage can generate new dam- 
age) , to an increase A K of the average connectivity, and 
to some function of the system size N. Actually, it can- 
not be directly proportional to N, because nodes that 
are part of the frozen core always remain undamaged 
asymptotically; hence, a rough upper limit is given by 
the number of nonfrozen nodes, which at K c scales as 
TV 2 / 3 . A lower bound can be derived by the number of 
relevant nodes, that almost certainly propagate damage, 
i.e. 7V 1//3 21 1 . To summarize, we approximate 



Ad ks c(N) N a d AK, 



with 1/3 < a < 2/3; replacing Ad and AK with differ- 
entials and integrating yields 



d(K, N) w cx (N) exp [c 2 (N) N a K] . 



(8) 



In simulations, we find a w 0.42, which is well within the 
range we expect from our theoretical considerations as 
discussed above. 

We now apply this dependence to obtain high-accuracy 
fits of Eq. ©in the interval 1.6 < K < 2.1 (Fig. H lower 
panels); these fits yield 

{Kf BN , df BN ) = (1.87 ± 0.04, 0.65 ± 0.05) (9) 

for RBN and, correspondingly, 

{Kf™, df TN ) = (1.725 ± 0.035, 0.52 ± 0.04) (10) 

for RTN. 

Interestingly, K s is close to, but distinct from the crit- 
ical connectivities Kf BN = 2 and K? TN = 1.845, as 
predicted by mean-field theory. However, a natural com- 
parison has to consider possible deviations of Kc at the 
SP limit from these values. An intuitive definition of crit- 
icality for finite N can be formulated in terms of marginal 
damage spreading. If at time t one bit is flipped, one re- 
quires at time t + 1 [7I 



d(t+l) = ( Ps )(K c )K c 



1, 



(11) 



where (p s )(K) is the average damage propagation prob- 
ability. Naturally, the iteration of this map implies d = 1 
for all t. Note that the relation: (p s )(K c )K c = 1 is exact 
only in the framework of the AA. In the SP limit, we 
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instead have to set the right hand side of Eq. [8] to unity; 
inversion then leads to 

xr-(N) « (12) 

Fig. |H shows K^ arse (N), using the values Ci(N),c 2 (N) 
obtained from numerical fits of Eq. [8] for both RBN and 
RTN. We find that both systems, in a very good approx- 
imation, obey the scaling relationship 

K s c parse (N) « b ■ N~ s + K™ (13) 

with b = 3.27±0.79, 5 = 0.85±0.07 and = 1.90375 ± 
0.005 for RBN and b = 3.853 ± 0.76, 6 = 0.736 ± 0.05 
and A c °° = 1.75598 ± 0.005 for RTN. Hence, in the limit 
N — » oo, we have 

k co,rbn = ! 90375 ± 0.005 (14) 
for RBN, and for RTN 

k oc,rtn = ! 75598 ± 0.005. (15) 

Thus for both RBN and RTN, in the SP limit and for 
N — > oo, the critical connectivity is considerably below 
the value K c predicted by the annealed approximation 
(notice that for small N < 128, however, K s c parse (N) > 

j£ annealed^ 

It is beyond the scope of this letter to discuss possible 
causes for these deviations in detail (this will be accom- 
plished in a long paper). In simulations, however, we 
find that the statistical distributions of damage sizes in 
the SP limit are highly skewed, with most configurations 
leading to vanishing damage, and a fat tail of large dam- 
age events. These skewed distributions imply that with 
finite sampling, we always underestimate d, and hence 
the true K^ parse (N) and K s will deviate even stronger 
from the A A. Also, local fluctuations in damage propa- 
gation cannot be neglected in this limit, as it is assumed 
in mean- field approaches. 

Discussion. We investigated finite size scaling of dam- 
age spreading in both RBN and RTN near the sparse 
percolation (SP) transition. We find that the average 
damage d, quantified in terms of the Hamming distance 
of initially nearby system states, scales cx over 
the whole range of sparse connectivities < K < 3 
studied in this letter. The scaling exponents 7 show a 
cross-over from negative to positive values at character- 
istic points K^ BN and Kf TN below the critical points 
j^RBN an( j j^RTN _ -y^r e estimated the critical connectiv- 
ities K^ parse (N) based on marginal damage spreading, 
and found deviations from the annealed approximation 
even in the limit of large N. Interestingly, recent stud- 
ies suggest that gene regulatory networks appear to be 
in the ordered regime and reside sli ght ly below the phase 
transition between order and chaos [151 ] , while theory had 
proposed the critical line to be an evolutionary attractor 



[HQ. Our study may contribute a possible explanation 
to these observations: in finite networks, scaling of dam- 
age with increasing N (e.g., as a consequence of gene 
duplications) can be expected to be under strict selec- 
tive control. Hence, the need for robust systems might 
drive evolution to K s rather than to K c . At the same 
time, however, K s is close enough to criticality to enable 
rich dynamical behavior, as required by biological cells. 
Certainly this question of where in the dynamical phase 
space biological networks reside cannot be answered con- 
clusively, given the current state of systems biology. Ex- 
perimental studies of regulatory networks and further 
studies of in-silico evolutionary processes that adapt dy- 
namical networks, e.g., RBNs or RTNs, to robustly solve 
specific computational and functional tasks are required. 
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